home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libconv / TRY / check2d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  18.2 KB  |  727 lines
  1. #include <stdio.h>
  2. #include <sys/time.h>
  3. #include <math.h>
  4. #include "conv.h"
  5.  
  6. /* --------- The following definitions change 
  7.         according to precision required -------- */
  8.  
  9.  
  10. #ifdef SINGLE       /* complex single precision    */
  11.  
  12. typedef float  this_variable;
  13. typedef float  this_real;
  14.  
  15. #define OPS_PER_ITER    2
  16.  
  17. void  sfir2d_(), simple_sfir2d_(), ornl_sfir2d_(), sfirm1d_(), scorm1d_();
  18. #define SIMPLE_FIR        simple_sfir2d_
  19. #define SIMPLE_IIR        simple_siir2d_
  20. #define SIMPLE_COR        simple_scor2d_
  21. #define ORNL_FIR        ornl_sfir2d_
  22. #define ORNL_IIR        ornl_siir2d_
  23. #define ORNL_COR        ornl_scor2d_
  24. #define MY_FIR            sfir2d_
  25. #define MY_IIR            siir2d_
  26. #define MY_COR            scor2d_
  27. #define MY_C_FIR        sfir2d
  28. #define MY_C_IIR        siir2d
  29. #define MY_C_COR        scor2d
  30. #define MY_M_FIR        sfirm1d_
  31. #define MY_M_IIR        siirm1d_
  32. #define MY_M_COR        scorm1d_
  33. #define MY_INIT            sinit_
  34. #define TOLERANCE        1.e-4
  35. #define THIS_IS_REAL
  36. #endif
  37. #ifdef DOUBLE            /* real double precision    */
  38.  
  39. typedef double  this_variable;
  40. typedef double  this_real;
  41.  
  42. #define OPS_PER_ITER    2
  43.  
  44. void  dfir2d_(), ornl_dfir2d_(), simple_dfir2d_(), dfirm1d_(), dcorm1d_();
  45. #define SIMPLE_FIR        simple_dfir2d_
  46. #define SIMPLE_IIR        simple_diir2d_
  47. #define SIMPLE_COR        simple_dcor2d_
  48. #define ORNL_FIR        ornl_dfir2d_
  49. #define ORNL_IIR        ornl_diir2d_
  50. #define ORNL_COR        ornl_dcor2d_
  51. #define MY_FIR            dfir2d_
  52. #define MY_IIR            diir2d_
  53. #define MY_COR            dcor2d_
  54. #define MY_C_FIR        dfir2d
  55. #define MY_C_IIR        diir2d
  56. #define MY_C_COR        dcor2d
  57. #define MY_M_FIR        dfirm1d_
  58. #define MY_M_IIR        diirm1d_
  59. #define MY_M_COR        dcorm1d_
  60. #define MY_INIT            dinit_
  61. #define TOLERANCE        1.e-10
  62. #define THIS_IS_REAL
  63.  
  64. #endif
  65. #ifdef COMPLEX            /* complex single precision    */
  66.  
  67. typedef complex this_variable;
  68. typedef float  this_real;
  69.  
  70. #define OPS_PER_ITER    8
  71.  
  72. void  cfir2d_(), ornl_cfir2d_(), simple_cfir2d_(), cfirm1d_(), ccorm1d_();
  73. #define SIMPLE_FIR        simple_cfir2d_
  74. #define SIMPLE_IIR        simple_ciir2d_
  75. #define SIMPLE_COR        simple_ccor2d_
  76. #define ORNL_FIR        ornl_cfir2d_
  77. #define ORNL_IIR        ornl_ciir2d_
  78. #define ORNL_COR        ornl_ccor2d_
  79. #define MY_FIR            cfir2d_
  80. #define MY_IIR            ciir2d_
  81. #define MY_COR            ccor2d_
  82. #define MY_C_FIR        cfir2d
  83. #define MY_C_IIR        ciir2d
  84. #define MY_C_COR        ccor2d
  85. #define MY_M_FIR        cfirm1d_
  86. #define MY_M_IIR        ciirm1d_
  87. #define MY_M_COR        ccorm1d_
  88. #define MY_INIT            cinit_
  89. #define TOLERANCE        1.e-4
  90. #define THIS_IS_COMPLEX
  91.  
  92. #endif
  93. #ifdef ZOMPLEX            /* complex double precision    */
  94.  
  95. typedef zomplex this_variable;
  96. typedef double  this_real;
  97.  
  98. #define OPS_PER_ITER    8
  99.  
  100. void  zfir2d_(), ornl_zfir2d_(), simple_zfir2d_(), zfirm1d_(), zcorm1d_();
  101. #define SIMPLE_FIR        simple_zfir2d_
  102. #define SIMPLE_IIR        simple_ziir2d_
  103. #define SIMPLE_COR        simple_zcor2d_
  104. #define ORNL_FIR        ornl_zfir2d_
  105. #define ORNL_IIR        ornl_ziir2d_
  106. #define ORNL_COR        ornl_zcor2d_
  107. #define MY_FIR            zfir2d_
  108. #define MY_IIR            ziir2d_
  109. #define MY_COR            zcor2d_
  110. #define MY_C_FIR        zfir2d
  111. #define MY_C_IIR        ziir2d
  112. #define MY_C_COR        zcor2d
  113. #define MY_M_FIR        zfirm1d_
  114. #define MY_M_IIR        ziirm1d_
  115. #define MY_M_COR        zcorm1d_
  116. #define MY_INIT            zinit_
  117. #define TOLERANCE        1.e-10
  118. #define THIS_IS_COMPLEX
  119.  
  120. #endif
  121.  
  122. this_real my_val[3] = {0., -.33, 1.};
  123.  
  124. /* ---------- The rest is the same    ---------------- */
  125.  
  126. #define MAX_SIZE        33
  127. #define DEFAULT_TRIALS        99
  128. #define MAX_INC            5
  129. #define MAX_TIMES        3
  130.  
  131.  
  132. #define ABS(a) ( ((a)>0) ? (a) : -(a))
  133. #define MAX(a,b) ( ((a) > (b)) ? (a) : (b))
  134.  
  135. void (*simple_function)();
  136. void (*ornl_function)();
  137. void (*my_function)();
  138. void GetArguments();
  139.  
  140. double check_this();
  141.  
  142. int parallel;
  143. int is_random;
  144. int n_trials;
  145. int len = 4;
  146.  
  147.     int size, n_trials, n_times, i_error;
  148.     int nf, ng, nh;
  149.     int nfx, ngx, nhx, nfy, ngy, nhy;
  150.     int ldf, nfx1, nfx2, nfy1, nfy2;
  151.     int ldg, ngx1, ngx2, ngy1, ngy2;
  152.     int ldh, nhx1, nhx2, nhy1, nhy2;
  153.     int incf, incg, inch;
  154.     this_variable *vx, *vy, *vz, *vt;
  155.     this_variable alpha, beta;
  156.  
  157.     double t, x, y, z, u;
  158.  
  159. main(argc,argv)
  160. int argc;
  161. char *argv[];
  162. {
  163.     int i, j, k;
  164.  
  165. /* ******************************************************* */
  166.     GetArguments( argc, argv);
  167. /* ******************************************************* */
  168.  
  169.     i_error = 0;
  170.  
  171.   for( k = 0; k < n_trials ; k++) { 
  172.  
  173.     get_values();
  174.  
  175. #ifdef THIS_IS_REAL
  176.     printf("%3d<>%1dx( (%3d,%3d)(%3d,%3d)(%3d,%3d) | (%3d,%3d)(%3d,%3d):(%3d,%3d)(%3d,%3d):(%3d,%3d)(%3d,%3d) | %5.2f, %5.2f ): ",
  177.         k,n_times, nfx,nfy, ngx,ngy, nhx,nhy,
  178.         nfx1,nfx2,nfy1,nfy2,
  179.         ngx1,ngx2,ngy1,ngy2,
  180.         nhx1,nhx2,nhy1,nhy2,
  181.         alpha, beta  );
  182. #endif
  183. #ifdef THIS_IS_COMPLEX
  184.     printf("%3d<>%3dx( (%3d,%3d)(%3d,%3d)(%3d,%3d) | (%3d,%3d)(%3d,%3d):(%3d,%3d)(%3d,%3d):(%3d,%3d)(%3d,%3d) ): ",
  185.         k,n_times, nfx,nfy, ngx,ngy, nhx,nhy,
  186.         nfx1,nfx2,nfy1,nfy2,
  187.         ngx1,ngx2,ngy1,ngy2,
  188.         nhx1,nhx2,nhy1,nhy2  );
  189. #endif
  190. /*
  191.     printf("%1dx(%3d,%3d,%3d | %2d,%3d,%3d | %2d,%3d,%3d | (%5.2f,%5.2f) | %2d,%3d,%3d | (%5.2f,%5.2f)): ",
  192.         n_times, nf,ng,nh, incf,nf1,nf2, incg,ng1,ng2, alpha.real,alpha.imag, inch,nh1,nh2, beta.real,beta.imag);
  193. */
  194.     fflush(stdout);
  195.  
  196.     vx = (this_variable *)my_malloc(((nf+1)) * sizeof( this_variable));
  197.     vy = (this_variable *)my_malloc(((ng+1)) * sizeof( this_variable));
  198.     vz = (this_variable *)my_malloc(((nh+1)) * sizeof( this_variable));
  199.     vt = (this_variable *)my_malloc(((nh+1)) * sizeof( this_variable));
  200.  
  201.     if( (vx == (this_variable *)0) || 
  202.     (vy == (this_variable *)0) ||
  203.     (vt == (this_variable *)0) ||
  204.     (vz == (this_variable *)0))  {
  205.     fprintf( stderr, "Malloc problem ... Exiting");
  206.     exit(0);
  207.     }
  208.     i = MY_INIT( &nf, vx);
  209.     j = MY_INIT( &ng, vy);
  210.  
  211.     do_it();
  212.  
  213.     printf("\n", x);
  214.  
  215.     my_free ( vx);
  216.     my_free ( vy);
  217.     my_free ( vz);
  218.     my_free ( vt);
  219.     if( i_error)
  220.     return (-i_error);
  221.   }
  222.   return (0);
  223. }
  224.  
  225.  
  226. do_it() 
  227. {
  228.     int ii, jj, k;
  229.     int ifx1, ifx2, ify1, ify2;
  230.     int igx1, igx2, igy1, igy2;
  231.     int ihx1, ihx2, ihy1, ihy2;
  232.     int this_ldf, this_ldg, this_ldh, inc;
  233.     this_variable *this_vx, *this_vy, *this_vz, *this_vt;
  234.     this_variable this_alpha, this_beta;
  235.  
  236. /*
  237. *   Checking against a simple convolution
  238. */
  239. #ifdef THIS_IS_REAL 
  240.     this_alpha = 1.0;
  241.     this_beta = 0.0;
  242.  
  243.         vy[0] = 1.0;
  244. #endif
  245. #ifdef THIS_IS_COMPLEX
  246.     this_alpha.real = 1.0;
  247.     this_alpha.imag = 0.0;
  248.     this_beta.real = 0.0;
  249.     this_beta.imag = 0.0;
  250.  
  251.         vy[0].real = 1.0;
  252.         vy[0].imag = 0.5;
  253. #endif
  254.     inc = 1;
  255.  
  256.     ifx1 = 1; ifx2 = nfx;
  257.     ify1 = 2; ify2 = nfy;
  258.     igx1 = 0; igx2 = ngx;
  259.     igy1 = 0; igy2 = ngy;
  260.     ihx1 = ifx1+igx1; ihx2 = nfx+ngx-1;
  261.     ihy1 = ify1+igy1; ihy2 = nfy+ngy-1;
  262.     this_ldf= ifx2;
  263.     this_ldg= igx2;
  264.     this_ldh= ihx2;
  265.     k = nh+1;
  266.     MY_INIT( &k, vz);
  267.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  268.     for( ii = 0; ii < n_times ; ii ++)
  269.     MY_FIR(        vx,&inc,&this_ldf,&ifx1,&ifx2,&ify1,&ify2,
  270.             vy,&inc,&this_ldg,&igx1,&igx2,&igy1,&igy2,
  271.             vt,&inc,&this_ldh,&ihx1,&ihx2,&ihy1,&ihy2,
  272.             &this_alpha,&this_beta);
  273.     for( ii = 0; ii < n_times ; ii ++)
  274.     SIMPLE_FIR( &ifx2, &ify2, &igx2, &igy2, vx, vy, vz);
  275.  
  276.     t = check_this( vz, vt, nh+1) ;
  277.  
  278. /*
  279. *   Checking IIR against a simple IIR_filter
  280. */
  281.     ifx1 = 1; ifx2 = nfx;
  282.     ify1 = 2; ify2 = nfy;
  283.     igx1 = 1; igx2 = ngx;
  284.     igy1 = 2; igy2 = ngy;
  285.     ihx1 = ifx1-igx1; ihx2 = nfx+ngx-1;
  286.     ihy1 = ify1-igy1; ihy2 = nfy+ngy-1;
  287.     k = nh+1;
  288.     MY_INIT( &k, vz);
  289.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  290.     for( ii = 0; ii < n_times ; ii ++)
  291.     MY_IIR(        vx,&inc,&this_ldf,&ifx1,&ifx2,&ify1,&ify2,
  292.             vy,&inc,&this_ldg,&igx1,&igx2,&igy1,&igy2,
  293.             vt,&inc,&this_ldf,&ihx1,&ifx2,&ihy1,&ify2);
  294.     for( ii = 0; ii < n_times ; ii ++)
  295.     SIMPLE_IIR( &ifx2, &ify2, &igx2, &igy2, vx, vy, vz);
  296.  
  297.     t = check_this( vz, vt, nh+1) ;
  298.  
  299. /*
  300. *   Checking COR against a simple corelation
  301. */
  302.     ifx1 = 0; ifx2 = nfx-1;
  303.     ify1 = 0; ify2 = nfy-1;
  304.     igx1 = 0; igx2 = ngx-1;
  305.     igy1 = 0; igy2 = ngy-1;
  306.     ihx1 = ifx1-igx2; ihx2 = nfx+ngx-1;
  307.     ihy1 = ify1-igy2; ihy2 = nfy+ngy-1;
  308.     this_ldf= nfx;
  309.     this_ldg= ngx;
  310.     this_ldh= nfx+ngx-1;
  311.     k = nh+1;
  312.     MY_INIT( &k, vz);
  313.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  314.     for( ii = 0; ii < n_times ; ii ++)
  315.     MY_COR(        vx,&inc,&this_ldf,&ifx1,&nfx,&ify1,&nfy,
  316.             vy,&inc,&this_ldg,&igx1,&ngx,&igy1,&ngy,
  317.             vt,&inc,&this_ldh,&ihx1,&ihx2,&ihy1,&ihy2);
  318.     for( ii = 0; ii < n_times ; ii ++)
  319.     SIMPLE_COR( &nfx, &nfy, &ngx, &ngy, vx, vy, vz);
  320.  
  321.     t = check_this( vz, vt, nh+1) ;
  322.  
  323. /*
  324. *   Checking the Fortran version against a Fortran template
  325. */
  326.     k = nh+1;
  327.     MY_INIT( &k, vz);
  328.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  329.     for( ii = 0; ii < n_times ; ii ++)
  330.     MY_FIR(        vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  331.             vy,&incg,&ldg,&ngx1,&ngx,&ngy1,&ngy,
  332.             vz,&inch,&ldh,&nhx1,&nhx,&nhy1,&nhy,
  333.             &alpha,&beta);
  334.     for( ii = 0; ii < n_times ; ii ++)
  335.     ORNL_FIR(   vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  336.             vy,&incg,&ldg,&ngx1,&ngx,&ngy1,&ngy,
  337.             vt,&inch,&ldh,&nhx1,&nhx,&nhy1,&nhy,
  338.             &alpha,&beta);
  339.     t = check_this( vz, vt, nh+1) ;
  340.  
  341.     k = nh+1;
  342.     MY_INIT( &k, vz);
  343.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  344.     for( ii = 0; ii < n_times ; ii ++)
  345.     MY_IIR(        vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  346.             vy,&incg,&ldg,&ngx1,&ngx,&ngy1,&ngy,
  347.             vz,&inch,&ldh,&nhx1,&nhx,&nhy1,&nhy);
  348.     for( ii = 0; ii < n_times ; ii ++)
  349.     ORNL_IIR(   vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  350.             vy,&incg,&ldg,&ngx1,&ngx,&ngy1,&ngy,
  351.             vt,&inch,&ldh,&nhx1,&nhx,&nhy1,&nhy);
  352.     t = check_this( vz, vt, nh+1) ;
  353.  
  354.     k = nh+1;
  355.     MY_INIT( &k, vz);
  356.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  357.     for( ii = 0; ii < n_times ; ii ++)
  358.     MY_COR(        vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  359.             vy,&incg,&ldg,&ngx1,&ngx,&ngy1,&ngy,
  360.             vz,&inch,&ldh,&nhx1,&nhx,&nhy1,&nhy);
  361.     for( ii = 0; ii < n_times ; ii ++)
  362.     ORNL_COR(   vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  363.             vy,&incg,&ldg,&ngx1,&ngx,&ngy1,&ngy,
  364.             vt,&inch,&ldh,&nhx1,&nhx,&nhy1,&nhy);
  365.     t = check_this( vz, vt, nh+1) ;
  366. /*
  367. *   Checking the C version against a Fortran version
  368. */
  369.     k = nh+1;
  370.     MY_INIT( &k, vz);
  371.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  372.     for( ii = 0; ii < n_times ; ii ++)
  373.     MY_C_FIR(   vx,incf,ldf,nfx1,nfx,nfy1,nfy,
  374.             vy,incg,ldg,ngx1,ngx,ngy1,ngy,
  375.             vt,inch,ldh,nhx1,nhx,nhy1,nhy,
  376. #ifdef THIS_IS_REAL
  377.             alpha,beta);
  378. #endif
  379. #ifdef THIS_IS_COMPLEX
  380.             &alpha,&beta);
  381. #endif
  382.     for( ii = 0; ii < n_times ; ii ++)
  383.     MY_FIR(        vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  384.             vy,&incg,&ldg,&ngx1,&ngx,&ngy1,&ngy,
  385.             vz,&inch,&ldh,&nhx1,&nhx,&nhy1,&nhy,
  386.             &alpha,&beta);
  387.     t = check_this( vz, vt, nh+1) ;
  388.  
  389.     k = nh+1;
  390.     MY_INIT( &k, vz);
  391.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  392.     for( ii = 0; ii < n_times ; ii ++)
  393.     MY_C_IIR(   vx,incf,ldf,nfx1,nfx,nfy1,nfy,
  394.             vy,incg,ldg,ngx1,ngx,ngy1,ngy,
  395.             vt,inch,ldh,nhx1,nhx,nhy1,nhy);
  396.     MY_IIR(        vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  397.             vy,&incg,&ldg,&ngx1,&ngx,&ngy1,&ngy,
  398.             vz,&inch,&ldh,&nhx1,&nhx,&nhy1,&nhy);
  399.     t = check_this( vz, vt, nh+1) ;
  400.  
  401.     k = nh+1;
  402.     MY_INIT( &k, vz);
  403.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  404.     for( ii = 0; ii < n_times ; ii ++)
  405.     MY_C_COR(   vx,incf,ldf,nfx1,nfx,nfy1,nfy,
  406.             vy,incg,ldg,ngx1,ngx,ngy1,ngy,
  407.             vt,inch,ldh,nhx1,nhx,nhy1,nhy);
  408.     for( ii = 0; ii < n_times ; ii ++)
  409.     MY_COR(        vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  410.             vy,&incg,&ldg,&ngx1,&ngx,&ngy1,&ngy,
  411.             vz,&inch,&ldh,&nhx1,&nhx,&nhy1,&nhy);
  412.     t = check_this( vz, vt, nh+1) ;
  413. /*
  414. *   Checking Mulitple 1D Filters along First dimension
  415. */
  416.     igy1 = 0;
  417.     igy2 = 1;
  418.     
  419.     k = nh+1;
  420.     MY_INIT( &k, vz);
  421.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  422.     for( ii = 0; ii < n_times ; ii ++)
  423.     MY_M_FIR(   vx,&incf,&ldf,&nfx1,&nfx,&nfy,
  424.             vy,&incg,&ngx1,&ngx,
  425.             vz,&inch,&ldh,&nhx1,&nhx,
  426.             &alpha,&beta);
  427.     for( ii = 0; ii < n_times ; ii ++)
  428.     MY_FIR(        vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  429.             vy,&incg,&ldg,&ngx1,&ngx,&igy1,&igy2,
  430.             vt,&inch,&ldh,&nhx1,&nhx,&nfy1,&nfy,
  431.             &alpha,&beta);
  432.     t = check_this( vz, vt, nh+1) ;
  433.  
  434.     k = nh+1;
  435.     MY_INIT( &k, vz);
  436.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  437.     for( ii = 0; ii < n_times ; ii ++)
  438.     MY_M_IIR(   vx,&incf,&ldf,&nfx1,&nfx,&nfy,
  439.             vy,&incg,     &ngx1,&ngx,
  440.             vz,&inch,&ldh,&nhx1,&nhx);
  441.     for( ii = 0; ii < n_times ; ii ++)
  442.     MY_IIR(        vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  443.             vy,&incg,&ldg,&ngx1,&ngx,&igy1,&igy2,
  444.             vt,&inch,&ldh,&nhx1,&nhx,&nfy1,&nfy);
  445.     t = check_this( vz, vt, nh+1) ;
  446.  
  447.     k = nh+1;
  448.     MY_INIT( &k, vz);
  449.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  450.     for( ii = 0; ii < n_times ; ii ++)
  451.     MY_M_COR(   vx,&incf,&ldf,&nfx1,&nfx,&nfy,
  452.             vy,&incg,     &ngx1,&ngx,
  453.             vz,&inch,&ldh,&nhx1,&nhx);
  454.     for( ii = 0; ii < n_times ; ii ++)
  455.     MY_COR(        vx,&incf,&ldf,&nfx1,&nfx,&nfy1,&nfy,
  456.             vy,&incg,&ldg,&ngx1,&ngx,&igy1,&igy2,
  457.             vt,&inch,&ldh,&nhx1,&nhx,&nfy1,&nfy);
  458.     t = check_this( vz, vt, nh+1) ;
  459. /*
  460. *   Checking Mulitple 1D Filters along Second dimension
  461. */
  462.     igx1 = 0;   
  463.     igx2 = 1;   
  464.  
  465.     k = nh+1;
  466.     MY_INIT( &k, vz);
  467.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  468.      for( ii = 0; ii < n_times ; ii ++)
  469.     MY_M_FIR(   vx,&ldf,&incf,&nfy1,&nfy,&nfx,
  470.             vy,&ldg,      &ngy1,&ngy,
  471.             vz,&ldh,&inch,&nhy1,&nhy,
  472.             &alpha,&beta);
  473.     for( ii = 0; ii < n_times ; ii ++)
  474.     MY_FIR(        vx,&incf,&ldf,&nfx1,&nfx ,&nfy1,&nfy,
  475.             vy,&incg,&ldg,&igx1,&igx2,&ngy1,&ngy,
  476.             vt,&inch,&ldh,&nfx1,&nfx ,&nhy1,&nhy,
  477.             &alpha,&beta);
  478.     t = check_this( vz, vt, nh+1) ;
  479.  
  480.     k = nh+1;
  481.     MY_INIT( &k, vz);
  482.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  483.      for( ii = 0; ii < n_times ; ii ++)
  484.     MY_M_IIR(   vx,&ldf,&incf,&nfy1,&nfy,&nfx,
  485.             vy,&ldg,      &ngy1,&ngy,
  486.             vz,&ldh,&inch,&nhy1,&nhy);
  487.     for( ii = 0; ii < n_times ; ii ++)
  488.     MY_IIR(        vx,&incf,&ldf,&nfx1,&nfx ,&nfy1,&nfy,
  489.             vy,&incg,&ldg,&igx1,&igx2,&ngy1,&ngy,
  490.             vt,&inch,&ldh,&nfx1,&nfx ,&nhy1,&nhy);
  491.     t = check_this( vz, vt, nh+1) ;
  492.  
  493.     k = nh+1;
  494.     MY_INIT( &k, vz);
  495.     bcopy( vz, vt, (nh+1)*sizeof(this_variable));
  496.      for( ii = 0; ii < n_times ; ii ++)
  497.     MY_M_COR(   vx,&ldf,&incf,&nfy1,&nfy,&nfx,
  498.             vy,&ldg,      &ngy1,&ngy,
  499.             vz,&ldh,&inch,&nhy1,&nhy);
  500.     for( ii = 0; ii < n_times ; ii ++)
  501.     MY_COR(        vx,&incf,&ldf,&nfx1,&nfx ,&nfy1,&nfy,
  502.             vy,&incg,&ldg,&igx1,&igx2,&ngy1,&ngy,
  503.             vt,&inch,&ldh,&nfx1,&nfx ,&nhy1,&nhy);
  504.     t = check_this( vz, vt, nh+1) ;
  505. }
  506.  
  507. void GetArguments( argc, argv)
  508. int argc;
  509. char *argv[];
  510. {
  511.     int i, j, k;
  512.     int nerror = 0;
  513.  
  514. #define ON    1
  515.  
  516.     parallel = 0;
  517.     is_random = 1;
  518.     n_trials = DEFAULT_TRIALS;
  519.  
  520. /* ******************************************************* */
  521.     for ( i = 1 ; (i < argc) && (nerror != ON) ; i ++ ) {
  522.     if( argv[i][0] == '-') {
  523.         switch ( argv[i][1]) {
  524.  
  525.         case 'i' :
  526.         case 'I' :  
  527.                 is_random = 0;
  528.                 break;
  529.         default  :  nerror = ON;
  530.         }
  531.     }
  532.     else {
  533.         if( i+1 > argc)
  534.         nerror = ON;
  535.         else { 
  536.         n_trials = atoi( argv[i]);
  537.         }
  538.     }
  539.     }
  540.     if( nerror == ON) {
  541.     fprintf( stderr, 
  542.         "Usage : %s [-p(arallel)] [n_trials]\n", argv[0]);
  543.     exit(0);
  544.     }
  545.     if( is_random)
  546.     srandom( (123*getpid()) | 0x01);
  547.     else 
  548.     srandom( 1);
  549. }
  550. get_values(){ 
  551.   int jj;
  552.     if( is_random) { 
  553.     nfx = random() % MAX_SIZE + 1;
  554.     ngx = random() % MAX_SIZE + 1;
  555.     nhx = nfx+ngx+random() % MAX_SIZE;
  556.  
  557.     nfx1 = random() % nfx;
  558.     nfx2 = nfx1 + random() % (nfx - nfx1);
  559.     ngx1 = random() % ngx;
  560.     ngx2 = ngx1 + random() % (ngx - ngx1);
  561.     nhx1 = nfx1+ngx1;
  562.     nhx2 = nhx1 + random() % (nhx - nhx1);
  563.  
  564.     jj = random() % 23 - 11;
  565.     nfx1 -= jj; nfx2 -= jj;
  566.     jj = random() % 23 - 11;
  567.     ngx1 -= jj; ngx2 -= jj;
  568.     jj = random() % 23 - 11;
  569.     nhx1 += jj; nhx2 += jj;
  570.     
  571.     nfy = random() % MAX_SIZE + 1;
  572.     ngy = random() % MAX_SIZE + 1;
  573.     nhy = nfy+ngy+random() % MAX_SIZE;
  574.  
  575.     nfy1 = random() % nfy;
  576.     nfy2 = nfy1 + random() % (nfy - nfy1);
  577.     ngy1 = random() % ngy;
  578.     ngy2 = ngy1 + random() % (ngy - ngy1);
  579.     nhy1 = nfy1+ngy1;
  580.     nhy2 = nhy1 + random() % (nhy - nhy1);
  581.  
  582.     jj = random() % 23 - 11;
  583.     nfy1 -= jj; nfy2 -= jj;
  584.     jj = random() % 23 - 11;
  585.     ngy1 -= jj; ngy2 -= jj;
  586.     jj = random() % 23 - 11;
  587.     nhy1 += jj; nhy2 += jj;
  588.     
  589.     incf = 1 + random() % MAX_INC;
  590.     incg = 1 + random() % MAX_INC;
  591.     inch = 1 + random() % MAX_INC;
  592.  
  593.     ldf = incf * nfx + random() % 3;
  594.     ldg = incg * ngx + random() % 3;
  595.     ldh = inch * nhx + random() % 3;
  596.  
  597.     nf = ldf * nfy;
  598.     ng = ldg * ngy;
  599.     nh = ldh * nhy;
  600.     
  601. #ifdef THIS_IS_REAL
  602.     alpha = my_val[random() % 3];
  603.     beta  = my_val[random() % 3];
  604. #endif
  605. #ifdef THIS_IS_COMPLEX
  606.     alpha.real = my_val[random() % 3];
  607.     alpha.imag = my_val[random() % 3];
  608.     beta.real  = my_val[random() % 3];
  609.     beta.imag  = my_val[random() % 3];
  610. #endif
  611.     n_times = random() % MAX_TIMES + 1;
  612.  
  613.     }
  614.     else {
  615.     printf( "Enter sizex of f : ");
  616.     scanf( "%d", &nfx);
  617.     printf( "Enter sizey of f : ");
  618.     scanf( "%d", &nfy);
  619.     printf( "Enter sizex of g : ");
  620.     scanf( "%d", &ngx);
  621.     printf( "Enter sizey of g : ");
  622.     scanf( "%d", &ngy);
  623.  
  624.     n_times = 1;
  625.  
  626.     nfx1 = 0;        nfx2= nfx-1;
  627.     nfy1 = 0;        nfy2= nfy-1;
  628.     ldf  = nfx;
  629.     nf   = ldf * nfy;
  630.  
  631.     ngx1 = 0;        ngx2= ngx-1;
  632.     ngy1 = 0;        ngy2= ngy-1;
  633.     ldg  = ngx;
  634.     ng   = ldg * ngy;
  635.  
  636.     printf( "Enter i0x of f : ");
  637.     scanf( "%d", &nfx1);
  638.     printf( "Enter i0y of f : ");
  639.     scanf( "%d", &nfy1);
  640.     printf( "Enter i0x of g : ");
  641.     scanf( "%d", &ngx1);
  642.     printf( "Enter i0y of g : ");
  643.     scanf( "%d", &ngy1);
  644.  
  645.     nhx1 = 0;    nhx2= nfx2+ngx2;
  646.     nhy1 = 0;    nhy2= nfy2+ngy2;
  647.     nhx = nhx2-nhx1+1;  nhy = nhy2-nhy1+1;
  648.     ldh  = nhx;
  649.     nh   = ldh * nhy;
  650.  
  651. #ifdef THIS_IS_REAL
  652.     alpha = my_val[1];
  653.     beta  = my_val[1];
  654. #endif
  655. #ifdef THIS_IS_COMPLEX
  656.     alpha.real = my_val[1];
  657.     alpha.imag = my_val[0];
  658.     beta.real  = my_val[0];
  659.     beta.imag  = my_val[1];
  660. #endif
  661.     incf = 1;
  662.     incg = 1;
  663.     inch = 1;
  664.     }
  665. }
  666.  
  667. double check_this( this_variable *a, this_variable *b, int n)
  668. {
  669.     double max, av;
  670.     int i, nerr=0;
  671.     max = av = 0.0;
  672.  
  673. #ifdef THIS_IS_REAL
  674.     for (i = 0; i < n; i++){
  675.     t = a[i] - b[i];
  676.     t = t * t;
  677.     if( t > max) max = t;
  678.     t = a[i];
  679.     av += t * t;
  680.     }
  681. #endif
  682. #ifdef THIS_IS_COMPLEX
  683.     for (i = 0; i < n; i++){
  684.     t = a[i].real - b[i].real;
  685.     z = a[i].imag - b[i].imag;
  686.     t = t * t + z * z;
  687.     if( t > max) max = t;
  688.     t = a[i].real;
  689.     z = a[i].imag;
  690.     av += t * t + z * z;
  691.     }
  692. #endif
  693.     max = sqrt(max);
  694.     av  = sqrt( av / (double)n);
  695.     if( av > 0.0) {
  696.     max = max / av;
  697.     }
  698.     if( max > TOLERANCE) {
  699.     fprintf(stderr, " Difference with reference is %.12e \n\n", max);
  700.     max = TOLERANCE * av;
  701.     max = max * max;
  702.     for (i = 0; i < n; i++){
  703. #ifdef THIS_IS_REAL
  704.         t = a[i] - b[i];
  705.         t = t * t;
  706. #endif
  707. #ifdef THIS_IS_COMPLEX
  708.         t = a[i].real - b[i].real;
  709.         z = a[i].imag - b[i].imag;
  710.         t = t * t + z * z;
  711. #endif
  712.         if( t > max) {
  713.         if( nerr < 10)
  714.             fprintf(stderr, " %d", i);
  715.         nerr++;
  716.         }
  717.     }
  718.     fprintf( stderr, "\n %d Errors detected \n", nerr);
  719.     }
  720.     else
  721.     printf("-", t);
  722.  
  723.      fflush(stdout);
  724.      i_error += nerr;
  725.     return(max);
  726. }
  727.